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This   paper  addresses  two   problems  of  interest   in   service 
system  analysis:     (a)   that  of  making  statistical,  data- 
driven   estimates  of  the  long-run  probability  of  a   long 
delay,  and  (b)   the  assessment  of  rate  of  approach  to  a 
long-run  system  performance  measure  such  as  expected 
delay,  the  rate  being  characterized  by  a   simple  exponen— * 
tial ,  at  least  initially.     Both  are  illustrated  by 
reference  to  M/G/l   and  related  systems. 

1 .  Introduction 

The  application  of  probability  theory  to  a  wide  variety  of  congestion  problems 
arising  in  communication  systems  has  been  well  catalogued  in  many  papers;  the 
treatises  by  Borovkov  (1984),  Cohen  (1969),  Cooper  (1972),  Cox  and  Smith  (1961), 
Feller  (1966),  Kleinrock  (1975),  Gross  and  Harris  (1974),  and  many  others  attest 
to  the  attractiveness  of  the  area  for  probability  modelling.   In  this  body  of  work, 
several  features  are  noticeable.  First,  most  of  the  elegant  solutions  obtained 
are  in  somewhat  implicit  form,  being  presented  as  functional  equations,  or,  fre- 
quently, as  integral  (Laplace)  transforms,  generating  functions,  and  sometimes  as 
combinations  of  the  above.  Moments  (particularly  the  first)  of  delays  or  number 
in  queue  under  long-run  conditions  when  steady  state  prevails  are  the  most  explicit 
and  assessible  performance  summaries.  Second,  the  results  obtained  are  presented 
in  terms  of  component  distribution  functions  and  stochastic  processes  (renewal, 
Poisson,  etc.)  that  are  taken  as  known;  only  rarely  are  issues  addressed  that 
arise  when  actual  data  is  to  be  used  as  a  basis  for  inference  from  the  models; 
however,  see  Cox  (1965).  Third,  only  fragmentary  comprehensible  information  con- 
cerning transient  behavior,  not  to  mention  time-dependence  of  parameters  is  avail- 
able; recently  however  work  by  Roth  and  Odoni  (1983),  Roth  (1981),  Abate  and 
Whitt  (1985),  Lee  (1985)  has  elucidated  the  simple  exponential -1  ike  approach  to 
the  steady  state  displayed  by  a  great  many  stochastic  systems  with  stationary, 
time-homogeneous  arrivals  and  service  processes. 

Instead  of  directly  dealing  with  the  above  classes  of  problems  much  attention 
has  been  concentrated  on  modelling  large  networks  of  servers,  particularly  steady- 
state  Markovian,  cf.  Kelly  (1979)  for  an  elegant  treatment;  see  package  programs 
created  by  AT&T  (the  D.  Mitra  PANACEA)  and  by  IBM  (work  by  Lavenberg,  et  al_) . 
Other  numerical  work  by  Neuts  (1981)  and  associates  has  pioneered  the  exploita- 
tion of  special  structure  of  certain  stochastic  systems.  And  a  considerable  ef- 
fort to  create  approximations,  ranging  from  heavy-traffic  and  diffusion  to,  lately, 
WKB  approaches,  see  Knessl  ,  Matkowsky  e_t  aj_  (1986)  has  been  evident;  see  in  par- 
ticular Newell  (1982)  for  pioneering  work. 

This  paper  deals  specifically  with  approaches  to  the  two  classes  of  somewhat 
neglected  problems  referred  to  earlier:  those  of  data-driven  inference,  and  of 
assessing  transient  behavior.  The  approaches  proposed  are  illustrated  concretely 

1 


in  terms  of  the  simple  M/G/l    system,   but  apply  more  widely. 

2.      Inference  Concerning  Long  Steady-State  Delays   in  M/G/1-Like  Systems 

Consider  a   single  service  system  approached  by  stationary  Poisson   (A)    traffic 
with  X   known.      Service  times,   or  message  lengths,    X,   are   independently  distributed 
as   Fx(x);   assume  AE[X]   =  p  <  1.     This   is  the  M/G/l    system.      Let  observations  of 
the  service  times   be  all    that  is   known  about   F:     denote  these  by  x-j  ,X2, . . .  ,xn. 
The  objective  is  to   supply  estimates  of  long-run  characteristics  of  the  system: 
in  particular,   estimates,   and  error  assessments   thereof,   of  the  probability  of  a 
long  delay  experienced  by  an  arriving  message. 

2.1     Virtual   Waiting  Time  Transform;  The  Long  Tail. 

It  is  well    known  that  if  W(t)   is  the  virtual   waiting  time  in  the  M/G/l   and 
p  <  1,  then  the  Laplace  transform 

E[e"sW]   =  Urn  E[e'sW(t)]   =  Np     gY         =  -^-  (2.1) 

t-  l-o[]-E[g"     J   1       ]-pA(s)  ~J 

1    PL     sE[X]        J 

where  A(s)   is  the  Laplace-Stiel tjes  transform  of  a  distribution.     It  can  also  be 
seen  that  in  various   interrupted  and  priority  service  systems,   including  those  in 
which  the  server  "takes  vacations,"  that  the  above  can  be  changed  to 

E[e'sW]     =         ]:p      .C(s)  /  (2.2) 

l-pB(s) 

where  possibly  B(s)  is  the  transform  of  a  completion  time  (effective  service  time, 
i.e.  X  including  durations  of  a  random  number  of  interruptions  occurring  therein), 
and  C(s)  is  the  transform  of  an  honest  distribution  that  accounts  for  effects  of 
busy  period  starts;  see  Gaver  (1968). 

If  A(s)  (or  B(s))  exists  for  s  >  s0,  s0  <  0  then  there  will  be  a  smallest 
real  zero,  -s  =  k  >  0,  of  the  denominator  of  (2.1 ),( (2.2) ) ,  which  can  be  used  to 
show  that 

P{W  >w}  -  D(<)e"KW,  w  -*•  °o  (2.3) 

One  way  of  establishing  the  above  exponential  tail  property  is  to  introduce 

^{s)=hm 1=  /  p{W>w}e"swdw 

s  0 

into   (2.1 )  which  leads  to 

Hs)   ^pC1^-]  +  pA(sH(s)  (2.4) 

equivalent  to  the  terminating  renewal    equation,   see  Feller   (1966),   p.   362, 

_  _  w 

Fu(w)   e  P{W>w}   =  H(w)    +  p     /     F(w-x)A(dx)  (2.5) 

0 

Introduce  ?w(w)  =  Fw(w)e6w,  6  real  and  positive,  into  (2.5) 

F?J(w)  =  H#(w)  +  /   F#(w-x)pe0XA(dx)  ,  (2.6) 

0 

and  choose  0  =   <  so  that 

00  oo 

/     pe9xA(dx)    =     /     A#(dx)    =  1    ,  (2.7) 

0  0 


— # 
yielding  a   standard  renewal    equation   for  F   ,   to  which  the   key  renewal    theorem 

applies,  yielding  as  w  ■*  °° 


•o 


/  H#(w)dw 
lim   F*(w)    =  lim   Fw(w)eKW  =   D(K)    =  5_ ,  (2.8) 

J  xA   (dx) 
0 

equivalent  to  (2.3).  A  similar  expression  is  available  for  (2.2). 

2.3  Statistical  Estimation  of  Probability  of  Long  Wait. 

If  data  are  available  on  service  times  (message  lengths)  then  an  estimate  of 
the  requisite  transform  is 

y    i  n   -sx. 
E[e  SX]  -\     I     e  '  (2.9) 

i=l 

and  k  is  the  unique  positive  solution  of  this  sample  equivalent  of  (2.7): 

n   ex 

I 
i=l 

A  three-term  Taylors'  series  expansion  around  6=0  gives  an  initial  estimate 


h\C-    I     e  i-l)   1  .  (2.10) 


k   =  2(-^-)  (2.11) 


7? 

T    k        k 
where  x  =  (x,  +  ...  +  x  )/n,  the  kth_  sample  moment;  (2.11)  is  recognized  as  being 

the  non-parametric  sample  version  of  the  exponential  distribution  parameter  ob- 
tained by  normalizing  to  W/E[W]  and  letting  p  t  1  (the  heavy  traffic  approximation). 
A  search  or  Newton-Raphson  iteration  starting  from  (2.11)  quickly  yields  <  as 
accurately  as  is  necessary. 

Some  theoretical  asymptotic  results  concerning  modes  of  behavior  of  k  will  now 
be  sketched;  more  detail  will  be  presented  elsewhere.  First,  it  is  essential  that 

E[eK  ]  <  *  in  order  for  (2.7)  to  hold,  and  so  by  continuity  f   (k)  =  E[X  eK  ]  <°° 
for  any  integer  m.  Now  express  (2.10)  as 

,   ,     n      9x. 
f(-)  =  Mtt    I    e       -1)1-  1  (2-12) 

9  n  1=1 

and  expand  in  Taylor's   series  around  <:     since  f(<)   =  0,   for  some  <£ 

0    =    f(K)     +    (K-K)f'(K)     +    ^-(K-K)2f"((|,K)      ,  (2,13) 


SO 

f(<) 


K  -  K 


f'(ic)     +    [(K-K)f"((f>K) 


(2,14) 


^ote  that  if  E[e2<X]   <  °°  then  Var[eKX]   <  °°  and 

2 
Var[f(K)]   ■  Vn  Var[eKX]    ,  (2.15) 

K 

ih  i  1  e 


E[f(K)]  =  0  . 

Additionally,  f  (k)  ■+  E[f  (k)]  by  the  weak  law  of  large  numbers  and  <-<  ■»■  0  in 
probability,  so  by  the  central  limit  theorem, 


/n(K-K)-  EEf'(^ -  N(0,1)  (2.16) 


2kX 
the  standard  Normal /Gaussian  distribution.  If,  however,  E[e   ]  =  °°,  then  if  any 

kind  of  limiting  distribution  is  to  exist,  the  components 

e 

of  the  sum  f(<)  must  be  in  the  domain  of  attraction  of  a  stable  law,  which  will  be 
assumed.  In  any  case  f^(<)  -»■  E[?'(k)]  as  before,  so  a  proper  normalization  by  a 
power  of  n  shows  that  (k-k)  has  a  stable  law  as  limiting  distribution.  Fortunate- 
ly, the  relatively  well-behaved  Normal/Gauss  limit  prevails  when  traffic  intensity 
is  relatively  high  (p  >  1/2  in  case  X  ~  Exp(u),  for  example).   ~.   % 

2.4  Assessing  Variability  of  the  Estimates. 

Having  obtained  an  estimate  of  the  parameter  k,  and  of  the  probability  of  a 
waiting  time  exceeding  large  t,  which  involves  k,  it  is  desirable  to  appraise  the 
errors  involved.  Two  error  sources  are:  the  systematic  error  '(bias)  resulting 
from  fitting  an  incorrect  model,  and  the  effect  of  finite  sample  size  (random 
error).  The  latter  is  the  easiest  to  evaluate,  provided  the  underlying  model  is 
nearly  correct.  The  bias  issue  is  more  difficult  to  deal  with;  one  approach  is 
to  begin  by  fitting  an  elaborate,  multiparameter  model  and  then  to  check  for  the 
contribution  of  extra  parameters  in  a  prediction  context.  Since  this  paper  takes 
a  non-parametric  approach  to  estimation,  the  bias  resulting  from  an  incorrect 
specification  of  the  service  time  distribution  is  presumably  absent,  although  the 
assumptions  of  stationary  Poisson  arrivals,  and  independence  are  still  reflected 
In  the  basic  transform  (2.1),  and  hence  in  the  estimates.  Unfortunately,  classi- 
cal procedures  that,  for  example,  avail  themselves  of  the  asymptotic  approxima- 
tions of  properties  of  maximum  likelihood  parameter  estimates  to  estimate  sampling 
errors  are  no  longer  available  in  the  present  environment.  We  have  chosen  in- 
stead to  investigate  the  performance  of  several  modern  procedures  often  recom- 
mended for  obtaining  standard  errors  of  estimate  and  confidence  Intervals:  the 
jackknife  (Quenouille,  Tukey,  Miller,  Hinkley,  and  others)  and  the  bootstrap 
(Efron  and  associates  and  others).  These  methods,  particularly  the  bootstrap, 
which  involves  simulation,  are  computer  intensive,  but  are  the  only  alternatives 
currently  known  for  dealing  with  complex  situations  such  as  that  at  hand. 

The  Jackknife 

The  above  name  refers  to  a  procedure  originally  introduced  by  Quenouille  for 
bias  reduction  (1956),  and  adapted  by  Tukey  (1958)  to  obtain  approximate  confidence 
intervals.  Suppose  interest  is,  on  a  parameter  G  (e.g.  our  k,  or  some  function 
thereof)  that  is  estimated  by  0,  using  a  complex  calculation  from  data 
(x-j  ,X2,A. .  ,xn) ,  just  as  k  is.  The  idea  is  that  of  assessing  variability  by  recom- 
puting k  after  removing  independent  subgroups  of  data  of  equal  size,  and  then 
using  the  recomputed  k  values  to  estimate  a  variance,  which  is  in  turn  applied  to 
state  a  standard  error  or  a  two-sided  confidence  interval  that  contains  the  true 
6  with  specified  confidence.  A  few  details  follow;  for  more,  see  Efron  (1982) 
and  his  more  recent  work,  or  Mosteller  and  Tukey  (1977).  The  actual  calculation 
involves  splitting  n  into  g  disjoint  groups  of  size  m;  n  =  mg.  Then  calculate 
ehj)'  J  =  l>2,...,g:  the  estimate  of  G  that  omits  the  jth_  group.  Now  Tukey 
(but  not  Efron)  computes  pseudo-values 


yd  =ge-  (g-Ue^j, 


which  are  then  treated  as  independent:  use  y  =  6jK  or  the  point  estimate  of  9, 
and  its  approximate  variance  as 

s^/g  =  J  (yj-y)2/(g-D(g)  ■    I  (^.jj-^.jj^Cg-D/g* 

as   it  turns  out.     Tukey  recommends   referring  y  to  Student's   t  with   g-1    degrees  of 
freedom  to  obtain  confidence  limits.     As   the  name   "jackknife"   is   intended  to 
imply,   the  tool    is   inexact  and  a   bit  crude  for  small    samples--just  as   a   true 
jackknife  is   not  well-adapted  for  delicate  surgery.      It   is   a   handy  non-parametric 
option. 

The  Bootstrap 

In  simplest  form  the  bootstrap  suggests  creating   from  observations   (x., 
i  =  l,...,n)   the  empirical   distribution 

a,         i     n  (  1 ,     x.    <  x 

Fn(x)   .1  J     l(Xl.x)s         I(x1tx).J  '"        . 

1=1  [  0,     x.    >  x 

1  S         _J 

Then  one  re-samples  wi£h  replacement  from  Fn  to  obtain  a  bootstrap  sample  ofsize 
n,   from  which  6   (e.g.   k(1))    is   calculated.      This^process   is   repeated  B   (e.g. 
300-500)   times,  and  the  empirical    distribution,    F~,  of  the   resulting  estimates, 
{<(b),  b  =  1,2,..,B}  is  employed  as  an  approximation  to  the  sampling  distribution 
of  k:     the  upper  and  lower  10%  points   of  F~  approximate  two-sided  80%  confidence 
limits   for  k,   for  instance.  / 

3.     Transient   Response:     Exponential    Characterizations 

Applications  of  queueing  theory  frequently  require  assessment  of  service  sys- 
tem transient  behavior,   particularly  when  the  system   is   initially  idle  and  heavy 
traffic  conditions  prevail,   for  then  approach  to  steady-state  conditions  may  be 
slow.     Odoni   and   Roth   (1983)    have  reviewed  work  on   the  nature  of  the  approach   to 
steady-state  conditions,   and  have  investigated  characterizations  of  that  approach 
that  are  simply  stated  and  comprehensible,   being  exponential .      In  this   section 
we  elaborate  upon  the  theme  of  exponential   approach.     Work  by  Abate  and  Whitt 
(1985)   has  similar  aims,  but  proceeds  somewhat  differently. 

3.1      From  Transform  to  Exponential 

Suppose  (X(t),   t>0}   is  the  state  variable  of  a  process,  and  that  expressions 


for 


OO  00 

•St^r  .,  .,,  ,  .  v    Iw,„xt   ,.  r        "St 


(a)     >Us)   =     /  e"il-E[^(X(t))|X(0)]dt  =     /e-iV(t)dt 
A  0  0  A 


(3.1) 


(b)  lim  E|>(X(t)  |X(0)]   <  °o 

are  known.     An  approximate  representation  of 

(c)  ty(t)   =  E[ip(X(t))|X(0)] 

1s  desired.  The  proposed  representation  is  a  simple  exponential  interpolation  of 
the  form 

Jx(t)  =  ^x(O)e~0t  +  u,x(oo)(i  -e"6t)  ,  (3.2) 

where  3  governs   the  rate  of  approach  to   the  steady-state  value.     The   proposal    is 
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to  identify  3  by  minimizing  the  integrated  squared  error 

00 

MB)  =  /  tyx(t)  -^x(0)e"3t+^x(c»)(l  -e"et)}2w(t)dt  .  (3.3) 

The  weight  function  w(t)  may  be  made  to  focus  upon  values  of  g   appropriate  for 
t-regions  of  interest.     To  optimize  on  3>  study 

CO 

^-  ■  K     /  tyx(t)   -ipx(0)e-3t  -  ^x(co)(l  -e"et)}te"etw(t)dt  (3.4) 

=  0   . 
This   is  equivalent  to 


oo 


/  .Ut)w(t)te"3tdt  =  iUO)   j  e'23ttw(t)dt  +  iU°°)   j  (l-e~et)e~6ttw(t)dt         (3.5) 
0       a  a      0  x      0 

If  w(t)  =  1  the  above  can  be  expressed  in  terms  of  the  derivative  of  jthe  Laplace 
transform  $»(s) : 

-  ^*x(8)  4tX(0)  t  |  ♦,(-)]  ,  (3.6) 

P 

often  a  transcendental  equation  that  must  be  solved  numerically  for  3.  Clearly, 
if  no  solution  or  multiple  solutions  exist  then  the  simple  universal  form  (3.2) 
is  called  into  question,  but  in  such  cases  weight  functions  are  usefully  invoked. 
If  interest  centers  on  ultimate  approach  (t  ■*  «)  to  ^x(°°)  then  the  smallest  3 
satisfying  (3.6)  is  of  interest. 

3.2  The  M/G/l  Example. 

As  an  important  illustration  of  the  above  ideas,  let  X(t)  ■  W(t)  be  the 
virtual  waiting  time  in  an  M/G/l  system,  and  consider  ijty(t)  =  E[W(t)|W(0)  =0]. 
We  know  that 

\,M   -^Poo'5'-?  (3'7) 

s 

where 

Poo(s)  =  s+A[l-b(s)]  =I1|;[l-b(s)]  =?h  (3-8) 

b(s)  being  the  transform  of  a  busy  period  duration,  satisfying 

b(s)  =  Fx(s  +  A(l-b(s))  ; 

Fx  is  the  Laplace-Stiel tjes  transform  of  the  service  time  or  message  length  dis- 
tribution Fx.-    Numerical   solution  of  (3.6)    is  now  possible.     An  approximation  to 
the  (smallest)   value  of  3  governing  the  ultimate  approach  (t  ->  °°)  when  p  t  1,   i.e. 
in  the  heavy-traffic  situation,   is 


V4     '"p       -h'"(0) 

In  general,   h'"(0)   will    involve  higher  moments  of  the  service  time;   in  case 
X  -  Exp(y)   considerable  simplification  occurs   and 

8  *  y(l-p)2  {8(l+p)}"0'5  (3 J0) 


i/hich  will  not,  of  course  agree  exactly  with  formulas  or  numerical  results  given 
jreviously  by  other  authors,  e.g.  see  Odoni  and  Roth  (1983)  being  a  compromise 
["Procrustean")  interpolation  of  a  single  exponential  over  the  infinite  t-range. 
[n  order  to  focus  upon  a  8-value  relevant  to  ultimate  approach  it  has  been  found 
^hat  a  weight  procedure  of  the  form  w,+,(t)  =  [l-exp(-B(k)t)]  is  workable  and 
iractable;  start  with  3(0)  =  0  (the  unweighted  procedure  described),  then  intro- 
iuce  w.+,  to  determine  B.+,  and  iterate  to  convergence,  which  occurs  rapidly. 

\  similar  approach  should  apply  to  find  a  6  appropriate  for  small  t. 

The  procedures  described  above  should  be  applicable  in  the  data-driven  context 
)f  the  problem  of  Section  1:  a  first  step  is  to  introduce  the  empirical  transfo 


rm 


:x  and  from  it  proceed  to  b(s),  to  ^u(s),  and  then  to  8'  Progress  in  this  area, 
ind  details  of  the  above  investigations  will 


be  reported  elsewhere 
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